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ABSTRACT 

This  paper  describes  an  implementation  of  broadband  soil 
attenuation  in  finite -difference  time  domain  (FDTD)  simulations 
of  seismic  wave  propagation  from  impulsive  sources.  We 
concentrate  on  the  attenuation  phenomenon,  the  computational 
approach,  comparisons  with  results  from  non-attenuating  soil 
models,  and  the  impacts  on  range  estimation.  Results  are  based 
on  our  three  dimensional  viscoelastic  FDTD  code  ptop,  which 
allows  impulsive  and  moving  vehicle  simulations  over  realistic 
heterogeneous  geologies  and  surface  topographies. 

Soil  attenuation  refers  to  the  decay  of  seismic  energy  by  intrinsic 
material  losses  in  soil.  It  can  reduce  the  amplitude  of  propagating 
waves  and  shift  the  frequency  of  signal  energy,  thus  affecting 
vehicle  range  estimates  and  seismic  signatures.  Relative  to  other 
attenuation  factors,  i.e.,  geometric  spreading  and  scattering,  soil 
attenuation  often  dominates  the  total  attenuation  characteristics 
of  seismic  surface  waves  generated  by  moving  vehicles. 

ptop  represents  soil  attenuation  using  the  viscoelastic 
constitutive  equations  of  a  “standard  linear  solid”  (SES).  The 
SES  form  implemented  is  a  material  symbolized  by  a  spring  and 
dashpot  in  series  that  are  in  parallel  with  a  spring.  The  viscous 
behavior  of  a  standard  linear  solid  can  be  characterized  across  a 
limited  frequency  band  by  the  quality  factors  Qp  and  Qs,  which 
are  material  properties  that  are  inversely  proportional  to  the 
attenuation  coefficients  of  P-  and  S-waves,  respectively. 

For  the  frequency  range  of  interest  to  seismic  vehicle  tracking, 
soils  are  known  to  have  an  approximately  constant  g.  To  model 
broadband  attenuation  in  our  computations,  each  Q  is  defined  in 
ptop  by  superimposing  the  band  limited  damping  mechanisms  of 
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L  standard  linear  solids  at  L  selected  frequencies  so  that  the 
corresponding  Q  is  approximately  constant  over  a  desired 
bandwidth.  This  results  in  a  g  function  that  is  consistent  with 
observed  attenuation  in  seismic  surface  waves. 

We  demonstrate  ptop^s  viscoelastic  soil  attenuation  response 
through  simulations  of  an  impulsive  source  at  the  surface  of  a 
realistic  geology  and  show  how  the  source’s  spectral  energy  is 
shifted  and  absorbed.  We  further  comment  on  the  impact  that 
soil  attenuation  has  on  target  range  estimates. 


1.  INTRODUCTION 


The  US  Army  is  interested  in  the  characteristics  of  seismic  surface  waves  from  moving  armored 
vehicles.  Systems  such  as  Hornet,  Raptor,  and  Rattler,  and  unattended  ground  sensor  networks  of 
the  Future  Combat  System,  will  rely  in  part  on  measured  seismic  surface  waves  for  successful 
operation.  One  feature  of  these  waves  is  that  they  show  smooth  amplitude  decay  as  a  function  of 
vehicle  range  (Moran  and  Greenfield,  1997,  and  Prado,  1998).  As  a  consequence,  tracking 
algorithms  for  estimating  range -to-target  using  measured  seismic  signals  are  generating 
increasing  interest  (Moran  et  ah,  1998). 

In  support  of  the  Army’s  seismic  sensing  needs  and  its  development  of  tracking  algorithms,  we 
are  developing  a  high-fidelity  numerical  model  for  seismic  wave  propagation  that  considers  the 
complex  effects  of  topographical  features,  shallow  geological  structure,  and  material  attenuation 
on  propagating  waves.  Its  results  can  be  used,  for  example,  in  lieu  of  field  data  for  developing, 
testing  and  refining  tracking  algorithms  for  networked  seismic  sensors  (Moran  et  ah,  2001). 

Soil  attenuation  refers  to  the  decay  of  seismic  energy  in  propagating  waves  by  intrinsic  material 
losses  in  soil.  It  can  reduce  the  amplitude  of  propagating  waves  and  shift  the  frequency  of  signal 
energy,  thus  affecting  vehicle  range  estimates  and  seismic  signatures.  Attenuation  in  near-surface 
soil  is  the  focus  of  this  paper.  Specifically  we  describe  soil  attenuation  as  a  phenomenon;  our 
approach  for  including  soil  attenuation  in  our  seismic  wave  propagation  model  and  in  the 
geologic  models  required  as  input;  a  demonstration  of  numerical  seismic  propagation  with  soil 
attenuation;  and  the  impacts  of  soil  attenuation  on  range  estimation. 

2.  SOIL  ATTENUATION 


The  propagation  of  small-strain  seismic  waves  in  soil  is  a  small  disturbance  phenomenon  that 
does  not  alter  the  fabric  of  the  soil  or  cause  permanent  deformations  (Santamarina  et  ah,  2001). 
Nonetheless,  energy  of  the  small  disturbances  can  be  dissipated  in  heat,  resulting  in  considerable 
changes  in  seismic  waveforms.  This  dissipation  is  the  loss  phenomenon  we  call  soil  attenuation. 

Soil  attenuation  is  included  in  one  of  the  three  classes  of  attenuation  mechanisms  in  seismic  wave 
propagation:  material,  geometric,  and  apparent  attenuation.  Material  attenuation  refers  to  the 
dissipative  loss  phenomena  of  geologic  materials,  including  soil.  It  occurs  almost  entirely  in  shear 
by  small  transverse  movements  of  material  lattices  and  grain  boundaries  (Lay  and  Wallace, 

1995).  Geometric  attenuation  is  the  spreading  of  wave  energy  as  it  propagates  outward  from  a 
source.  For  surface  waves,  geometric  spreading  is  cylindrical.  Apparent  attenuation  includes 
scattering  of  energy  by  heterogeneous  geologic  and  topographic  features  and  partial  energy 


transmission  across  geologic  interfaces  (Santamarina  et  al.,  2001).  Apparent  attenuation  is  the 
most  difficult  to  quantify.  Relative  to  the  other  classes,  material  attenuation — i.e.,  soil 
attenuation — often  dominates  the  total  attenuation  characteristics  of  seismic  surface  waves 
generated  by  moving  vehicles. 


For  a  Rayleigh  surface  wave,  the  dependence  of  the  particle  velocity  amplitude.  A,  on  range,  r, 
from  the  source,  is  of  the  form  (Lay  and  Wallace,  1995) 


exp 


7^ 

Qc 


(1) 


where  Ao  is  a  reference  particle  velocity  amplitude;  g  is  a  unitless  “quality  factor”  that  quantifies 
the  material  attenuation  and,  depending  on  the  data,  the  apparent  attenuation  as  well;  and / and  c 
are  the  frequency  and  wave  speed  of  the  propagating  wave,  respectively.  As  Equation  1  indicates, 
lower  Q  results  in  higher  attenuation.  Constant  Q,  over  an  applicable  frequency  bandwidth,  is 
recognized  to  be  the  most  physically  realistic  functional  variation  of  Q  to  represent  material 
attenuation  in  seismic  wave  propagation  (e.g.,  Xu  and  McMechan,  1998). 

Using  Equation  1,  Greenfield  and  Moran  (2001)  processed  field-measured  data  sets  over  limited 
frequency  bands  to  perform  least-squares  particle  velocity  amplitude  versus  range  fits  of  the  data. 
These  fits — from  five  Army  proving  ground  sites — are  shown  as  amplitude  in  dB  vs.  range 
curves  in  Figure  1.  The  strong  influence  of  geologic  structure  from  site  to  site  is  evident.  The 
figure  shows  large  differences  in  both  amplitude  levels  and  in  attenuation  rate.  The  different 
amplitude  levels  are  due  principally  to  differences  in  material  stiffness  at  the  sites.  The  different 
attenuation  rates,  where  there  is  overlap  of  the  ranges  of  the  data  in  the  curves,  can  be  due  to 
differences  in  both  material  and  scattering  losses.  Based  upon  the  geology  of  the  sites,  however, 
differences  in  soil  losses  are  likely  to  have  had  far  greater  influence  in  the  slopes  of  the  Figure  1 
curves  than  differences  in  scattering  losses.  Thus,  as  the  Figure  1  curves  suggest,  synthetic  data  to 
be  used  by  researchers  for  developing  range  estimation  algorithms  should  include  the  effects  of 
soil  attenuation  on  propagating  waves. 


3.  COMPUTATIONAL  APPROACH 


Our  work  (e.g.,  Moran  et  al.,  1999,  Ketcham  et  al.,  2000)  modeling  seismic  wave  propagation 
originated  from  the  linear  elastic  formulations  and  numerical  implementations  described  by 
Flestholm  and  Ruud  (1998),  who  incorporated  surface  topography  with  an  appropriate  stress-free 
surface  boundary  condition  into  a  finite-difference  time  domain  (FDTD)  wave  propagation  model 
featuring  8*’'-order,  staggered-grid,  finite -difference  operators.  To  accommodate  surface 
topography,  Flestholm  and  Ruud  express  geologic  models  using  a  curvilinear  grid  that  is 
transformed  into  a  rectangular  computational  grid  of  equal  grid  spacing.  This  mapping  can  be 
visualized  by  proportionally  stretching  the  rectangular  grid  in  the  vertical  direction  so  that  the  free 
surface  matches  the  topographic  function.  Appendix  A  provides  a  partial  mathematical 
description  of  the  Flestholm  and  Ruud  (1998)  elastic  formulation.  Specifically  it  describes  the 
curved  system,  equations  of  motion,  and  stress-strain-velocity  relationships. 

To  model  soil  absorption  in  finite  difference  seismic  wave  propagation  calculations,  geophysical 
researchers  commonly  use  the  theoretical  framework  of  linear  viscoelasticity.  Relative  to  linear- 
elastic  materials,  Xmsw-viscoelastic  materials  retain  the  linearity  between  stress  and  strain,  but  the 
relationship  also  includes  time.  A  linear-elastic  body  has  a  simple  memory — it  “remembers”  only 
its  unstrained  configuration,  allowing  the  present  deformation  to  be  found  knowing  only  the  load. 
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Figure  1.  Vertical  particle  velocity  amplitude  vs.  range  fits  of  data  from  Army  proving  ground  sites. 


Figure  2.  A  “standard  linear  solid.”  /t'is  viscosity,  [ij  and  are  elastic  moduli,  and  cris  stress. 
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Figure  3.  Relaxation  and  creep  loading  solutions  to  the  stress-strain  relation  of  the  standard  linear 
solid.  The  solid  has  material  properties  of  the  upper  soil  layer  in  the  Figure  4  geologic  model. 


For  a  linear-viscoelastic  body,  however,  the  deformation  cannot  be  determined  unless  the  entire 
history  of  loading  is  known  (Fung,  1965). 

A  linear-viscoelastic  material  is  one  for  which  infinitesimal  strains,  displacements,  and  velocities 
are  considered  and  the  stress  tensor  (7y(x,y,z,t)  is  related  to  the  strain  tensor  eki{x,y,z,t)  by  a 
convolution  integral.  This  may  be  abbreviated  in  the  equivalent  forms 


tf//  G  ijki  dSki/  dt  Ski  dGijki/  dt 


(2) 


where  Gyu  is  the  tensorial  relaxation  function  of  the  material  and  *  denotes  convolution  integral 
operations  (Fung,  1965).  For  an  isotropic  material,  Gya  can  be  defined  by  two  scalar  relaxation 
functions,  one  in  isotropic  compression.  A,  and  one  in  shear,  2  M,  such  that 


(Jy  dJ\-/  dt  ^  3y  Skk 


2  dM/  dt  *  Sy  . 


(3) 


Because  the  relaxation  functions  are  scalar  functions,  it  is  appropriate  to  specify  their  behavior 
using  analogies  of  common  networks  of  springs  and  dashpots.  Our  FDTD  seismic  wave 
propagation  code,  ptop,  implements  viscoelasticity  according  to  the  finite  difference  formulations 
of  Flestholm  (1999)  and  Robertsson  et  al.  (1994).  These  formulations  use  the  constitutive 
equations  of  a  “standard  linear  solid”  (SLS  )  in  the  form  symbolized  by  a  “Maxwell  element” 

(i.e.,  a  spring  and  dashpot  in  series)  that  is  in  parallel  with  a  spring.  This  is  illustrated  in  Figure  2. 
In  this  figure  is  the  viscosity  of  the  solid,  i.e.,  the  proportionality  constant  between  stress  and 
the  rate  of  strain,  and  pi  and  jl2  are  elastic  moduli  of  the  solid,  i.e.,  proportionality  constants 
between  stress  and  strain.  p2  is  the  relaxed  modulus  corresponding  to  the  long-time  viscous 
response,  while  pi  +  p2  is  the  immediate,  unrelaxed  modulus.  The  stress-strain  relation  is  (e.g., 
Fung,  1965) 


(7+  Tedo/dt  =  P2  (£+  ds  /  dt  ) 


(4) 


where  and  Tg  are  the  relaxation  times  for  stress  and  strain,  respectively,  which  are  defined  by 
Ta  =  p'iPi  +  P2)/piP2,and  (5) 

Te  =  p'/pi.  (6) 

Solutions  to  the  stress-strain  relation  are  illustrated  in  Figure  3  for  a  Poisson  solid,  i.e.,  a  material 
with  Poisson  ratio  of  0.25.  The  solutions  are  for  relaxation  and  creep  loading.  They  have  the 
forms 


(7(0  =  [  1  -  (  1  -  Te/Ta)  exp(  -t /Za)  ]  1(0  , 

for  relaxation  loading — i.e.,  a  unit  step  application  of  strain,  and 
(7  (0  =  10^2  [  1  -  (  1  -  Ta/Tf  )  exp(  -t  /Ze)  ]  1(0  , 


(7) 


(8) 


for  creep  loading — i.e.,  a  unit  step  application  of  stress.  1(0  is  the  unit  step  function.  Shear  and 
compressive  responses  are  depicted  in  Figure  3  for  corresponding  shear  and  compressive  unit 
disturbances  that  are  shown  at  the  bottom  of  the  figure.  These  disturbances  are  applied  and 
removed  as  indicated.  The  responses  shown  are  for  a  viscoelastic  soil  with  shear  wave  speed  = 
577  m/s  and  density  =  1750  kg/m^,  one  material  considered  in  a  geologic  model  to  be  described 
subsequently. 

With  reference  to  Equation  3,  Robertsson  et  al.  (1994)  and  Flestholm  ( 1 999)  define  FI  =  A  -i-  2  M, 
and  use  the  standard-linear-solid  relaxation  solutions 


n  (0  =  (A  +  2  ^)  [  1  -  (1  -  zf/Za)  exp(  -t/Za )  ]  1(0  ,  and 
M(t)  =  p[  \  -  Ze/Za)  exp(  -t/z„)  ]  1(0  , 


(9) 

(10) 


to  define  material  behavior  in  their  finite  difference  formulations.  Here  A  and  n  are  Lame’s  elastic 
parameters,  t/  and  t/  are  the  strain  relaxation  times  for  P-  and  S-waves,  respectively,  and  is 
the  stress  relaxation  time  for  both  P-  and  S-waves.  Using  these  equations  Robertsson  et  al.  (1994) 
and  Hestholm  (1999)  derived  the  viscoelastic  stress-strain- velocity  relations  that  are  presented  in 
Appendix  B.  These  relations,  relative  to  the  elastic  stress-strain-velocity  relations  in  Appendix  A, 
include  the  effects  of  relaxation  times  and  memory  variables  on  the  stress.  The  memory  variables, 
Kij,  relate  the  current  value  of  stress  to  strain  history,  which  is  required  for  viscoelastic  modeling. 
As  indicated  in  Appendix  B,  Equations  9  and  1 0  can  be  generalized  to  represent  the  relaxation 
behavior  of  superpositioned  SLS  mechanisms. 

From  Figure  3  one  can  note  that  the  relaxation  of  a  Poisson  standard  linear  solid  is  dominated  by 
shear  behavior,  which  is  consistent  with  soil  attenuation  phenomena.  Indeed,  SFS  spring  and 
dashpot  constants  can  be  assigned  so  that  soil  absorption  during  seismic  wave  propagation  can  be 
realistically  modeled.  The  SFS  spring  and  dashpot  configuration  is  commonly  used  because  the 
viscous  behavior  can  be  characterized  across  a  limited  frequency  band  by  the  quality  factors  Qp 
and  Qs,  which  are  inversely  proportional  to  the  attenuation  of  P-  and  S-waves,  respectively. 

Qp  and  Qs  represent  the  attenuation  of  P-  and  S-waves  in  the  following  way;  they  refer  to  the 
fractional  loss  of  energy  per  cycle  of  a  P-  or  S-wave  oscillation  (Fay  and  Wallace,  1995), 
respectively.  The  inverse  proportionality  to  the  attenuation  coefficients  cCp  or  a,  is  given  by  \/Qp 
=  UpXp/n  for  P-waves  and  l/Q^  =  a,Xs/7rfor  S-waves,  where  Xp  and  A,  are  the  respective 
wavelengths.  Qp  and  Qs  can  also  be  defined  (Carcione  et  al.,  1988,  Xu  and  McMechan,  1998)  in 
the  frequency  domain  using  the  real  and  imaginary  parts  of  the  unrelaxed  bulk  and  shear  moduli. 
For  the  unrelaxed  shear  modulus,  jXu,  (which,  at  the  instant  of  loading  is  analogous  to  jd]  +  ld2  for 
the  SFS  in  Figure  2)  this  takes  the  MQ  form 

a“‘(/-)  =  lm[^„(/-)]/Re[At„(/-)],  (11) 

where /is  frequency,  jdu,  for  a  single  SFS  mechanism,  is  defined  as  (Xu  and  McMechan,  1998) 

Ai„(/-)  =  2^  [(  1  +  f  2;r/T/)  /  (  1  +  /  2;r/T/)] ,  (12) 

where  jx  is  the  relaxed  (elastic)  shear  modulus.  Defining  1  /InZa  as  the  relaxation  frequency/; , 
g/'  for  a  single  SFS  mechanism  is  then  (Xu  and  McMechan,  1998) 


Xu  and  McMechan  (1998)  similarly  define  Qp~^  and  the  unrelaxed  bulk  modulus.  Their  focus, 
however,  is  not  on  single  SFS  mechanisms  but  on  the  superposition  of  multiple  SFS  mechanisms 
in  order  to  achieve  nearly-constant  Q  to  realistically  represent  material  attenuation  over  the 
dominant  frequency  bandwidth  of  a  given  FDTD  calculation.  For  L  mechanisms.  Equations  12 
and  13  become  (Xu  and  McMechan,  1998) 
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Again,  Xu  and  McMechan  (1998)  similarly  define  Qp~^  and  the  unrelaxed  bulk  modulus.  They 
further  describe  an  inversion  technique  to  calculate  appropriate  relaxation  times  t/  and  t/  for 
each  mechanism.  Their  technique,  automated  in  software  (Ramos-Martinez  et  ah,  2000),  uses 
simulated  annealing  to  find  the  combination  of  L  relaxation  mechanisms  to  best  fit  the  desired  Q. 
The  relaxation  frequencies  are  chosen  according  to  the  desired  bandwidth  where  Q  is  to  be  nearly 
constant  and  are  assigned  identically  for  both  P-  and  S-wave  propagation. 

Because  the  unrelaxed  moduli  are  functions  of  frequency,  the  P-  and  S-wave  speeds  of  the 
material  vary  with  the  frequency  of  oscillating  disturbances.  Typically,  near  surface  geologic 
materials  get  stiffer  with  increasing  frequency  of  loading.  As  one  would  expect  from  elastic 
material  property  definitions,  P-  and  S-wave  speeds  increase  with  frequency  as  well.  An 
expression  for  the  unrelaxed  shear  wave  speed,  Vsu,  with  frequency  is  (Carcione  et  al.,  1988) 

r.„(/-)  =  Re[(At„(/-)/p)'^"].  (16) 

A  similar  expression  can  be  found  for  the  unrelaxed  compression  wave  speed. 

Stress  calculations  using  viscoelastic  relations  in  Appendix  B  rather  that  the  elastic  relations  in 
Appendix  A  are  the  principal  difference  between  viscoelastic  and  elastic  results  from  ptop.  As 
commented  in  Appendix  B,  the  calculations  can  be  performed  for  L  SLS  mechanisms.  For  this 
case,  an  approach  to  define  a  material’s  viscoelastic  properties  can  be  outlined  by: 

•  Select  the  desired  constant  Qp  and  Qs  to  define  soil  attenuation. 

•  Select  relaxation  frequencies of  L  standard  linear  solid  mechanisms  to  span  the 
bandwidth  where  constant  Q  is  desired. 

•  Find  the  P-wave  and  S-wave  strain  relaxation  times  t/  and  t/  of  each  mechanism  that 
provide  near-constant  values  of  Q,  approximating  constant  Qp  and  Qs,  over  the  selected 
bandwidth. 

This  is  the  approach  used  for  the  viscoelastic  material  definitions  reported  in  this  paper.  We 
adopted  the  Xu  and  McMechan  (1998)  simulated  annealing  technique  (Ramos-Martinez  et  al., 
2000)  to  find  the  strain  relaxation  times  of  each  mechanism. 

4.  DEMONSTRATION  OF  PTOFS  VISCOELASTIC  SOIL 

ATTENUATION 


Finite-difference  seismic  simulations  over  the  expected  range  of  coverage  for  FCS  seismic 
unattended  ground  sensor  systems  require  substantial  models  and  computational  durations.  As  a 
consequence,  our  simulations  are  parallel  computations  based  upon  a  domain  decomposition 
strategy  (Moran  et  al.,  1999).  We  perform  computations  on  multi-processor  computers  available 
at  DoD  High  Performance  Computing  Modernization  Program  centers.  The  simulations  described 
in  this  paper  were  performed  using  a  Cray  T3E  at  the  Army  High  Performance  Computing 
Research  Center.  The  geologic  model,  which  is  illustrated  in  Figure  4,  has  typical  topographic 
and  geologic  features  of  the  “hill  and  dale”  class  of  research  and  development  models  we  use  for 


testing  the  numerical  robustness  of ptop.  Its  size  is  roughly  350  m  by  500  m  by  80  m  (deep).  A 
2.8-m  node  spacing  defines  the  rectangular  grid  spacing  in  the  W-E  and  S-N  directions.  A 
minimum  of  1.6  m  is  the  node  spacing  in  the  vertical  direction,  which  is  elongated  to  handle 
topography  as  defined  in  Appendix  A.  The  side  and  bottom  boundaries  were  set  with  a  21 -cell- 
thick  absorption  layer  (Cerjan  et  ah,  1985),  providing  a  reduction  of  at  least  -20  dB  in  reflecting 
wave  particle  energy.  The  simulations  were  tests  of  the  long-duration  stability  of  our  viscoelastic 
code.  The  simulated  durations  were  up  to  12  s  with  time  steps  of  180  p,s. 

4. 1  Geologic  Model  and  Impulsive  Source 

The  geologic  model  is  a  synthetic  model  consisting  of  two  fairly  stiff  soil  layers  (above  and 
below  a  water  table)  overlying  granitic  bedrock.  It  is  similar  to  a  model  previously  reported  for 
demonstration  of  our  moving  impulsive  source  capability  (Ketcham  et  ah,  2000),  yet  larger  in 
horizontal  extent.  Two  common  geological  features  distinguish  its  sloping  topography:  an 
outcropping  of  the  bedrock  and  a  trench  representative  of  an  eroded  streambed.  Figure  4a  is  a 
surface  contour  graph  illustrating  the  topography  and  these  features. 

The  outcrop  is  roughly  elliptical  with  dimensions  of  80  m  by  200  m.  Its  peak  is  offset  laterally 
from  the  center  of  the  streambed  by  approximately  200  m.  The  streambed  is  roughly  100  m  wide 
by  8  m  deep.  “Downhill”  on  the  model  is  from  North  to  South,  i.e.,  from  the  top  of  Figure  1  to  the 
bottom,  as  a  0.002  slope  occurs  over  the  model  in  this  direction. 

Figure  4b  illustrates  the  subsurface  layering  of  the  model;  it  is  a  slice  at  the  185-m  South-North 
coordinate.  The  shades  of  the  model  refer  to  different  materials.  The  upper  two  layers  away  from 
the  outcrop  are  the  soil  layers.  The  surface  soil  layer  is  approximately  11m  thick.  The  lower  soil 
layer — i.e.,  the  soil  beneath  the  “water  table” — is  approximately  15  m  thick.  The  actual  values 
vary  throughout  the  model,  as  the  surface  is  not  flat.  In  addition,  both  soil  layers  reduce  in 
thickness  adjacent  to  the  outcrop  due  to  the  increasing  elevation  of  the  bedrock  surface  as  it  rises 
toward  the  outcrop,  and  the  upper  soil  layer  thins  toward  the  streambed. 

The  outcrop  features  an  upper  weathered  zone.  This  zone  is  depicted  in  Figure  4a  by  the  shading 
changes  beneath  the  outcrop.  The  uppermost  nodes  in  this  zone  have  seismic  propagation 
properties  identical  to  the  surface  soil  layer.  The  properties  vary  linearly  until  the  granitic  layer 
properties  are  reached. 

Table  1  lists  the  seismic  properties  of  the  three  principal  layers.  The  elastic  properties  P-wave 
speed,  S-wave  speed,  and  density  are  supplemented  by  the  quality  factors  Qp  and  Qs. 


Table  1.  Seismic  properties  of  layer  materials  in  geologic  model. 


Layer 

Compression- 
wave  speed  (m/s) 

Shear-wave 
speed  (m/s) 

Density 

(kg/m^) 

Material  loss  quality 
factors  Qp,  Q, 

Upper  soil 

1000 

577 

1750 

25,9 

Lower  soil 

1600 

625 

2000 

30,  15 

Granitic  bedrock 

3500 

2333 

2650 

150, 67 

Force  input  to  the  models  was  a  near-surface  vertical  pulse  with  a  maximum  force  of  1  N  and  a 
duration  of  approximately  0.03  s.  It  was  located  near  the  intersection  of  the  W-E  and  N-S  centers 
of  the  model,  just  west  of  the  trench,  at  W-E=252  m,  S-N=l  84.8  m,  and  elevation=70.3  m.  The 
depth  below  the  surface  was  7.2  m. 


0  50  100  150  200  250  300  350  400  450  500 

W-E[m] 


Figure  4.  (a)  Surface  contour  graph  of  geologic  model  with  weathered  bedrock  outcrop,  eroded 
streambed,  and  gently  sloping  flats,  (b)  Slice  through  model  at  South-North  coordinate  =  185  m. 
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Figure  5.  Spectra  of  Qp~^,  Q~^,  and  unrelaxed  P-  and  S-wave  speeds  for  the  1-mechanism  analysis 
tpld  fzl.  Solid,  dashdot,  and  dashed  curves  are  respectively  for  the  top-to-bottom  layers  in  Table  1. 
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Figure  6.  Spectra  of  and  unrelaxed  P-  and  S-wave  speeds  for  the  2-mechanism  analysis 

tplc  fzl.  Solid,  dashdot,  and  dashed  curves  are  respectively  for  the  top-to-bottom  layers  in  Table  1. 
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Figure  7.  Spectra  of  Q~'^,  Q~^,  and  unrelaxed  P-  and  S-wave  speeds  for  the  3-mechanism  analysis 
tplb  fzl.  Solid,  dashdot,  and  dashed  curves  are  respectively  for  the  top-to-bottom  layers  in  Table  1. 


4.2  Simulations  and  Results 


Three  viscoelastic  analyses  and  one  elastic  analysis  are  reported  here  for  the  geologic  model  and 
impulsive  source  described  above.  Table  2  presents  the  viscoelastic  specifications  that  distinguish 
the  analyses. 


Table  2.  Elastic  and  viscoelastic  simulations. 


Simulation 

ID 

Analysis  type 

Number  of  SLS 
mechanisms 

Relaxation 
frequencies  (Hz) 

tple  fzl 

Elastic 

0 

N/A 

tpldfzl 

Viscoelastic 

1 

40 

tple  fzl 

Viscoelastic 

2 

20,  200 

tplb  fzl 

Viscoelastic 

3 

2,  20,  200 

The  three  viscoelastic  analyses,  tpld  fzl,  tplc  fzl  andtplb  fzl,  use  1,  2  and  3  standard  linear 
solid  mechanisms,  respectively,  to  represent  the  material  attenuation  behavior.  The  relaxation 
frequency  of  the  1  -mechanism  analysis  corresponds  to  the  center  of  dominant  energy  of  the 
impulsive  source.  The  frequencies  of  the  2-  and  3 -mechanism  analyses  span  a  widening 
bandwidth  over  which  the  model  response  was  expected  to  be  substantial. 

Figure  5  illustrates  spectra  of  Qp~^  and  Qs~^  for  the  single-mechanism  analysis.  Curves  are  shown 
for  the  three  main  layers  of  the  geologic  model.  These  were  calculated  using  the  Xu  and 
McMechan  (1998)  inversion  technique  (Ramos-Martinez  et  ah,  2000).  The  peaks  of  the  Qp~^  and 
gC*  curves  in  Figure  5  occur  exactly  at  the  40-Flz  relaxation  frequency.  Individual  markers  in  the 
plots  show  the  desired  gp“'  and  gC*  values  at  the  relaxation  frequency.  Comparisons  of  the 
markers  with  the  spectra  reveal  that  the  inversion  fits  at  their  peak  are  accurate  for  the  lower  soil 
and  bedrock,  but  have  noticeable  error  for  the  upper  soil  with  the  highest  1/g  values.  Considering 
inaccuracies  that  occur  in  the  estimation  of  Q  from  field  data,  the  error  is  not  excessive. 

The  lower  three  plots  of  each  side  in  Figure  5  show,  from  top  to  bottom,  the  spectra  of  the 
unrelaxed  P-  and  S-wave  speeds  of  the  upper  soil  layer,  the  lower  soil  layer,  and  the  bedrock, 
respectively,  that  are  predicted  by  Equation  1 6  and  its  P-wave  counterpart.  The  increases  in  wave 
speeds  with  frequency  are  clear. 

Figures  6  and  7  show  spectral  quantities  corresponding  to  those  plotted  in  Figure  5  for  the  2-  and 
3 -mechanism  viscoelastic  analyses.  In  these  figures,  however,  the  frequency  bands  are 
increasingly  wider.  As  for  the  1  -mechanism  analysis  there  is  error  observed  in  the  spectral 
peaks  relative  to  the  desired  values,  but  the  error  is  considerable  mostly  for  gp“'  in  the  upper  soil. 
The  lower  plots  in  these  figures  again  show  the  increases  in  unrelaxed  wave  speeds  of  the  three 
layers  with  frequency.  Flere  the  increases  are  sustained  over  the  wider  bandwidths,  and  are 
thereby  greater. 

Figure  8  depicts  images  of  the  vertical  particle  velocity,  w,  on  the  surface  of  the  geologic  model. 
The  results  are  from  the  2-mechanism  viscoelastic  analysis.  These  are  shown  at  three  times 
during  the  propagation  of  seismic  waves  from  the  impulsive  source.  A  scale  in  the  figure  gives 
the  correspondence  between  the  image  shade  and  the  velocity  amplitude  in  m/s.  The  center  of  the 
scale  is  0  m/s.  Lighter  shades  indicate  positive/upward  velocities  and  darker  shades  indicate 
negative/downward  velocities.  The  principal  waveforms  displayed  in  the  images  are  fundamental 
Rayleigh  surface  waves,  which  have  cylindrical  decay  (1/r'^^,  r=radius  from  source)  in  the 
absence  of  the  topography  and  geology  that  disturb  this  decay.  The  waveforms  reveal  the 
geology-induced  diffraction,  refraction  and  reflection  of  the  waves. 
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Figure  8.  Images  of  vertical  particle  velocity  h>  on  geologic  model  surface  at  0.1  s,  0.2  s,  and  0.3  s. 
Image  shade  corresponds  to  value  of  w  in  scale.  Contour  lines  are  at  2-m  intervals. 
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Figure  9.  Location  of  three  lines  of  “receivers”  on  the  surface  of  the  geologic  model.  The  lines  are 
labeled  West,  East,  and  North  to  indicate  their  direction  from  the  source  epicenter. 


The  Figure  8  images  were  constructed  from  the  output  of  the  2-mechanism  simulation  at  each 
finite-difference  grid  point  on  the  surface,  providing  a  spatial  resolution  of  2.8  m.  Figure  9  shows 
another  image,  this  from  the  elastic  simulation,  with  three  lines  of  “receivers”  on  the  surface  of 
the  model.  The  lines  are  labeled  West,  East,  and  North  to  indicate  their  direction  from  the  source 
epicenter.  To  compare  the  elastic  and  viscoelastic  results,  w  versus  time  signals  were  collected 
from  the  full  output  of  the  analyses  at  the  illustrated  locations.  As  examples  of  the  early  w 
responses  from  the  simulations.  Figure  10  shows  elastic  and  2-mechanism  viscoelastic  signals 
from  three  of  the  locations:  one  at  the  epicenter  (Figure  10a),  one  along  the  north  line  at  50-m 
range  (Figure  10b),  and  another  along  the  north  line  at  100-m  range  (Figure  10c). 

The  large  early  features  of  the  Figure  10b  and  c  waveforms  are  the  Rayleigh  surface  waves.  The 
attenuation  relative  to  the  elastic  waveforms  and  the  effect  of  the  increasing  wave  speed  with 
frequency  shown  in  Figure  6  are  evident.  The  Rayleigh  wave  speed  for  the  upper  soil  layer  of  the 
elastic  analysis  is  530  m/s  (Rayleigh  wave  speed  Vr  =  0.92  for  this  Poisson  solid).  For  a  100-m 
range  this  would  arrive  at  0. 19  s,  which  is  consistent  with  what  the  100-m  range  signal  in  Figure 
10a  shows.  The  Rayleigh  wave  arrival  for  the  2-mechanism  analysis  is  earlier,  suggesting  an 
unrelaxed  surface  wave  speed  of  around  600  m/s  and  therefore  a  shear  wave  speed  of  650  m/s. 
Indeed,  the  shear  wave  speed  for  the  upper  soil  layer  in  Figure  6  is  close  to  650  m/s  in  the 
frequency  bandwidth  of  interest.  The  later,  oscillating  features  of  the  Figure  10b  and  c  waveforms 
appear  to  be  resonant  vibrations  associated  with  guided  energy,  i.e.,  energy  that  is  mostly  trapped 
in  the  soil  layers. 

Figure  11a  and  b  present  magnitudes  of  cross  spectra  between  the  epicenter  signals  of  Figure  10a 
and  the  50-  and  100-m  range  signals  in  Figures  10b  and  10c.  The  elastic  results  are  shown  on  the 
left  in  the  figure,  and  the  2-mechanism  viscoelastic  results  are  shown  on  the  right.  Shifts  in  the 
dominant  spectral  peaks  toward  higher  frequencies  from  the  cross  spectra  of  the  elastic  signals  to 
the  viscoelastic  signals  are  evident,  as  are  reductions  in  spectral  amplitude  relative  to  the 
reference  signal  at  0-m  range.  It  is  clear  that  spectral  response  quantities  are  affected  by  modeled 
soil  attenuation. 


5.  IMPACTS  ON  RANGE  ESTIMATION 


Figure  12  illustrates  maximum  signal  amplitude  (zero-to-peak  vertical  particle  velocity  w)  from 
the  three  geophone  lines  depicted  in  Figure  9.  The  results  in  each  plot  are  from  the  four  FDTD 
analyses  detailed  in  Table  2.  The  figure  highlights  the  impacts  that  attenuation  mechanisms  have 
on  simulated  amplitude  versus  range  data.  In  each  plot  the  elastic  result  has  the  highest  vertical 
particle  velocity  amplitudes  over  the  range  shown.  The  curves  reduce  noticeably  in  amplitude  as 
1-,  2-,  and  3 -mechanisms  are  added,  although  the  reduction  appears  greatest  when  adding  one  and 
then  two  mechanisms.  This  result  substantiates  the  previous  statement  made  upon  examination  of 
the  Figure  1  field  data  curves,  i.e.,  that  synthetic  data  to  be  used  by  researchers  for  developing 
range  estimation  algorithms  should  include  the  effects  of  soil  attenuation  on  propagating  waves. 

Figure  12  also  reveals  the  considerable  impact  that  geologic  and  topographic  features  can  have  on 
particle  velocity  amplitude  versus  range  data.  The  West  line  shows  a  sudden  drop  in  amplitude  as 
the  soil  thickness  becomes  small  near  the  bedrock  outcrop.  Figure  4  shows  the  section  of  the 
subsurface  geology  beneath  this  line.  It  is  at  a  range  of  just  over  100  m  where  the  amplitude  drops 
steeply,  which  is  just  what  one  would  expect  by  looking  at  the  Figure  4  section.  The  East  line 
shows  no  dramatic  change  as  the  trench  is  encountered,  but  relative  to  the  West  and  North  lines,  it 
is  evident  that  the  trench  has  the  effect  of  reducing  vertical  particle  velocity  amplitudes  at  ranges 
of  approximately  100  m  and  less. 


cross  spectrum,  magnitude  [dB] 


Figure  10.  Elastic  and  2-mechanism  viscoelastic  signals  from  three  “North  line”  receiver  locations: 
(a)  the  epicenter,  (b)  50-m  range,  and  (c)  100-m  range  . 
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Figure  11.  Cross  spectra  magnitudes  between  the  epicenter  signals  of  Figure  10a  and  the  50-  and 
100-m  range  signals  in  Figures  10b  and  10c.  (a)  Elastic  and  (b)  2-mechanism  viscoelastic  results. 
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Figure  12.  Maximum  signal  amplitudes  from  the  West,  East,  and  North  receiver  lines.  Results  from 
the  elastic  analysis  are  the  top  lines  with  circle  symbols.  The  1-,  2-,  and  3-mechanism  viscoelastic 
results  plot  from  top  (triangle  symbol)  to  bottom  (plus  symbol)  below  the  elastic  result,  respectively. 


Figure  13.  Sets  of  true  and  measured  amplitude  versus  range  curves  (top  plots)  and  estimated  range 
versus  true  range  curves  (bottom  plots)  for  three  hypothetical  field  sites  with  soil  attenuation  rates 
0.2  dB/m  (left  plots)  ,  0.1  dB/m  (middle  plots),  and  0.05  dB/m  (right  plots). 


Figure  12  makes  clear  that  attenuation  due  to  both  soil  losses  and  scattering  are  important  to 
model  when  producing  synthetic  results  for  developing  ground  sensor  range  estimation 
algorithms.  Another  issue — with  direct  field  implications — is  how  different  attenuation  rates  at 
different  sites  can  affect  estimates  of  range.  For  example,  consider  two  sites  that  have  identical 
particle  velocity  amplitude  measurements  in  dB  at  some  arbitrary  range,  but  that  have  soils  with 
significantly  different  attenuation  rates.  For  the  case  where  site  “one”  has  a  higher  attenuation 
rate  than  site  “two,”  the  errors  in  range  estimation  will  be  lower  at  site  one  due  simply  to  the 
steeper  slope  of  the  attenuation  curve  at  site  one.  In  Figure  1 3  we  present  results  of  an  exercise 
with  three  sites  to  illustrate  this. 

Figure  1 3  shows  sets  of  curves  for  three  hypothetical  field  settings  that  are  distinguished  solely  by 
the  soil  attenuation  rate  a.  The  three  rates  are  0.2,  0.1,  and  0.05  dB/m,  which  represent  a  realistic 
range  of  attenuation  rates  that  can  be  observed  in  field  tests  in  different  geologies. 

An  applicable  form  for  particle  velocity  amplitude  A  versus  range  r,  derivable  from  Equation  1 ,  is 
A=Ao  exp(-ar),  where  a  is  the  attenuation  rate  in  dB/m  and  Jo  is  a  reference  amplitude.  For 
the  three  settings  we  define  the  error  in  estimating  the  signal  amplitude  to  have  a  standard 
deviation  of  5  dB.  Perfect  amplitude  versus  range  data  are  generated,  which  are  depicted  by  the 
“true”  data  in  the  upper  plots  of  Figure  13,  and  normally  distributed  random  error  at  the  5  dB 
level  is  added.  The  data  with  error  are  used  to  estimate  the  range  as  if  they  were  measured  data, 
and  they  are  so  designated  in  Figure  13a,  c,  and  e.  The  estimated  ranges,  based  upon  a  range 
estimation  algorithm,  are  shown  in  the  bottom  plots  of  Figure  13  relative  to  the  true  ranges,  rms 
errors  in  the  range  estimates  are  listed  in  the  bottom  plots.  These  data  are  repeated  in  Table  3. 

The  rms  errors  decrease  with  increasing  attenuation  rate.  This  is  expected  because  as  the 
attenuation  rate  increases,  a  constant  error  in  the  measured  amplitude  in  dB  leads  to  a  smaller 
change  in  the  estimated  range. 

Table  3.  Range  estimation  errors  for  different  attenuation  rates. 


a  [dB/m] 

rms  error  in  range  estimate  [m] 

0.05 

64.2 

0.1 

35.8 

0.2 

24.3 

The  effect  of  soil  attenuation  on  range  estimation  becomes  more  complicated  in  the  presence  of 
background  seismic  noise  and  complex  geology.  Nonetheless,  where  signals  are  above  detectable 
levels,  range  estimates  are  likely  to  be  more  accurate  in  higher  soil  attenuation  areas. 

6.  CONCLUSION 


The  viscoelastic  seismic  simulations  we  have  presented  demonstrate  the  effect  of  soil  attenuation 
on  seismic  waveforms.  The  spectral  energy  from  an  impulsive  source  was  increasingly  absorbed 
and  shifted  as  the  number  of  relaxation  mechanisms  and  bandwidth  of  the  attenuation  was 
increased.  Because  soil  attenuation  affects  quantities  of  interest  for  range  estimation  and 
signature-spectra  calculations,  it  is  clear  that  synthetic  data  to  be  used  by  researchers  for 
developing  range  estimation  and  signature  identification  algorithms  should  include  the  effects  of 
soil  attenuation  on  propagating  waves. 

The  approach  used  for  defining  the  broadband  soil  attenuation,  i.e.,  the  Xu  and  McMechan  (1998) 
simulated  annealing  technique,  proved  an  effective  method  for  calculating  relaxation  times 
associated  with  desired  quality  factors  and  thereby  defining  the  viscoelastic  properties  of  our 


geologic  models.  Used  in  the  context  of  our  viscoelastic  seismic  wave  propagation  code  ptop,  we 
observed  expected  effects  of  Qp,  Qs,  and  unrelaxed  P-  and  S-wave  speeds  on  surface  waveforms 
and  spectra. 

Finally,  in  a  range  estimating  exercise  using  hypothetical  amplitude  versus  range  data  with 
different  attenuation  rates,  we  addressed  the  impact  of  soil  attenuation  on  range -to-target 
estimates.  We  observed  that  range  estimates  are  likely  to  be  more  accurate  in  higher  soil 
attenuation  areas  where  signals  are  above  detectable  levels. 
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APPENDIX  A.  ELASTIC  FDTD  FORMULATION 
(HESTHOLM  AND  RUUD,  1998) 


Curved  system  x,  y,  z  with  topography  funetion  (x,  y),  and  reetangular  system 


x  =  ^,  y  =  K,  z  =  -  Zo(x,y); 
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Partial  derivatives  needed  in  transformed  medium  equations  : 
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Equations  of  motion  in  {^,K,r\)  spaee : 
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Elastie  stress  -  strain  -  veloeity  relationships  in  {q,K,r])  spaee : 
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where :  p  is  density,  A  =  p(E^^  -  2V/)  and  p  =  pC/  are  Eame's  parameters  for  isotropic 
material  with  compression  wave  speed  and  shear  wave  speed  and 

are  the  components  of  the  body  force;  u,  v,  and  w  are  the  components  of  the  particle 
velocity;  and  a^^ ,  a^^ ,  a ,  a ,  a ,  and(7^^  are  the  stress  components. 


APPENDIX  B.  VISCOELASTIC  STRESS-STRAIN-VELOCITY 
RELATIONSHIPS  (HESTHOLM,  1999) 
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where  n  is  the  relaxation  modulus  for  P  -  waves  corresponding  to  A  +  2^7  in  the  elastic 
case  and  M  is  the  relaxation  modulus  for  S  -  waves  corresponding  to  jl  in  the  elastic  case; 
rf  andTg  are  the  strain  relaxation  times  for  P  -  and  S  -  waves,  respectively ;  is  the  stress 
relaxation  time  for  both  P  -  and  S  -  waves;  and  and  are 

memory  variables  relating  current  value  of  stress  to  strain  history.  For 
superposition  of  L  mechanisms,  summations  are  performed  in  the  equations 
above,  with  the  relaxation  moduli  defined  by 
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where  l(t)  is  the  Heaviside  or  unit  step  function. 


